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Abstract 

The generalized stacking fault (GSF) energy surfaces have received con- 
siderable attention due to their close relation to the mechanical properties 
of solids. We present a detailed study of the GSF energy surfaces of silicon 
within the framework of density functional theory. We have calculated the 
GSF energy surfaces for the shuffle and glide set of the (111) plane, and that 
of the (100) plane of silicon, paying particular attention to the effects of the 
relaxation of atomic coordinates. Based on the calculated GSF energy sur- 
faces and the Peierls-Nabarro model, we obtain estimates for the dislocation 
profiles, core energies, Peierls energies, and the corresponding stresses for 
various planar dislocations of silicon. 



1 Introduction 



An issue of central importance in materials science is the intrinsic ductility or 
brittleness of solids. To address this question, one has to consider how a crack 
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in the solid will respond to external loading. Since dislocation nucleation at a 
crack tip will cause the tip to blunt, it is customary to associate this process with 
ductile behavior. In contrast, brittle behavior is associated with crack propagation 
without dislocation emission, corresponding to the creation of a sharp crack tip. 
Therefore, the intrinsic ductility or brittleness of a solid can be determined by 
comparing the likelihood for dislocation nucleation at the crack tip to Griffith's 
criterion for cleavage (Griffith (1920)). 

Estimating the likelihood of dislocation nucleation at a crack tip under external 
loading is a non-trivial task. Significant advances have been made by Rice and 
collaborators recently in modehng dislocation nucleation at a crack tip based on 
the Pcicrls stress concept (Rice (1992), Bcltz and Rice (1991, 1992), Rice and Beltz 
(1994), Rice, Beltz, and Sun (1992), Beltz (1992), Sun, Beltz, and Rice (1993)) In 
that work, an important solid state parameter, the unstable stacking fault energy 
denoted by 7^5, was identified as the controlling parameter for dislocation emission 
at a crack tip under shear loading. One of the attractive features of Rice's theory 
is that the value of can be obtained from theoretical calculations without 
any ambiguity: Consider the process by which an infinite crystal is cut in half 
along a plane, and the upper half is sheared with respect to the lower half by 
a displacement vector f and let $(/) be the energy per unit area on the slip 
plane associated with this displacement. The energy surface $(/) obtained as a 
function of the generalized displacement vector / is called the generalized stacking 
fault (GSF) energy surface. The GSF energy surface or 7 surface and its relation 
to mechanical properties of solids were first considered by Vitek (1966), (1967), 
(1968), Vitek and Yamaguchi (1973), and Yamaguchi and Vitek (1975). In this 
context, Rice's unstable stacking energy 7„s is the lowest energy barrier in the 7 
surface that has to be surmounted during the shearing process that takes a crystal 
from an ideal configuration to another equivalent one. 

In addition to its relevance to dislocation emission, the GSF energy surface can 
be used to obtain the restoring force due to the misfit of the lattice near the core 
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of a dislocation (see for example, Joos, Ren, and Duesbery (1994)). By incorporat- 
ing these restoring forces into the Peierls-Nabarro model (Peierls (1940), Nabarro 
(1947)), one can calculate dislocation profiles, core energies, Peierls energies, and 
Peierls stresses for planar dislocations on the corresponding crystal plane. The 
calculation of dislocation properties within the Peierls-Nabarro model represents a 
drastic simplification, because it neglects the discreteness of the lattice by treating 
the problem in a continuum picture. Despite its obvious limitations, this approach 
provides a useful phenomenological framework for comparison of the structural 
properties as well as the energetics of various dislocations. 

From the above discussion, it is evident that an accurate description of the 
GSF energy surface is desirable. Empirical methods based on classical interatomic 
potentials are not sufficiently accurate for this task (Duesbery, Michel, Kaxiras, 
and Joos (1991)). Given the significance of GSF energy surfaces, it is important to 
obtain these values from a first-principles theoretical point of view, i.e. from calcu- 
lations which are free of adjustable parameters. First-principles calculations based 
on Density Functional Theory (DFT) (Hohenberg and Kohn (1964), Kohn and 
Sham (1965)), which are computationally much more demanding than empirical 
approaches, have been shown to be very successful in predicting the energy dif- 
ferences between different structures of solids. We have performed first-principles 
DFT calculations for both the shuffle and glide sets of the (111) plane as well as 
the (100) plane of silicon in the diamond lattice. There are several reasons for 
choosing these particular planes. These are low index planes of this crystal sys- 
tem, with in-plane dense packing of the atoms and relatively large spacing between 
planes. Therefore, these planes are natural candidates to be exposed during cleav- 
age. Moreover, the Burgers vectors associated with dislocations on these planes are 
short compared to other planes, which suggests a smaller dislocation core energy. 

The remaining of this paper is organized as follows: Section Q describes the 
computational techniques used in our first-principles calculations for GSF energy 
surfaces. Section | contains our results for the GSF energy surfaces of the (111) 
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and (100) planes of silicon. The effects due to the relaxation of atomic coordinates 
are also discussed there. In Section |^ we give a brief review of the Peierls-Nabarro 
model and then present the dislocation properties obtained within this model from 
the calculated GSF energy surfaces. We conclude with some remarks on the ap- 
plicability of Peierls-Nabarro model in Section ^. 

2 Computational Methods 

The local density approximation (LDA) to the exchange- correlation functional of 
DFT proposed by Perdew and Zunger (1984) was used for the GSF energy calcu- 
lations. The valence electron wave functions were expanded in a plane wave basis. 
The highest kinetic energy of the plane waves included in the basis set is 8 Ry. 
The ionic potential, including the screening by core electrons, was modeled by a 
nonlocal norm-conserving pseudopotential from Bachelet, Greenside, Barraff, and 
Schliiter (1982), and the scheme of Kleinman and Bylander (1982) was employed 
to make the potential separable in Fourier space. The d angular-momentum com- 
ponent was treated as the local part of the potential with the s and p components 
containing the nonlocal contributions. To simulate the block shearing process, a 
supercell consisting of 12 atomic planes in the direction perpendicular to the cut 
was used in these calculations. For the reciprocal space integration, we have used 
20 special k-points in the irreducible Brillouin zone, using the scheme of Monkhorst 
and Pack (1976). The approach of Car and Parrinello (1985), which allows simul- 
taneous relaxation of both the ionic and the electronic degrees of freedom, was 
employed in the present calculations. The minimum energy was obtained with the 
steepest descent method. Atoms on the four atomic planes farthest from the plane 
of the cut were kept frozen to simulate the effects of the bulk. Atomic relaxations 
are taken into account my minimizing the magnitude of forces calculated from the 
Hellmann-Feynman theorem. The structures are considered fully relaxed when the 
magnitude of the forces is smaller than 0.005 Ry/a.u. 
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The unrelaxed GSF energy surface is obtained by simply moving one half of 
the crystal rigidly with respect to the other half. This GSF energy surface can be 
altered when atoms on either side of the cut are allowed to respond to shearing 
forces. The difference in GSF energy surfaces with and without atomic relaxation 
can be very important, and can even change siginificantly the magnitude of as 
well as the configuration corresponding to it (sec next section). 

There is one subtlety that needs to be clarified in relation to the relaxation of 
ionic coordinates. Within Rice's continuum theory (Rice (1992)) the GSF energy 
is defined as a function of the displacement on either side of a mathematical cut 
in the middle of the crystal. In a real crystal, the closest approximation to this 
displacement is the relative displacement f of the two atomic planes immediately 
adjacent to the cut. When the two blocks on either side of the cut are moved rigidly, 

— * 

/ is identical to the relative displacement of the centers of the two halves of the 
crystal. When relaxation is included, special care must be taken in interpreting 

— * 

the proper value of / from the in-plane coordinates of the ions. For certain values 
of /, where $(/) is required by symmetry to have an extremum (as in the case of 
7us), only relaxations that are perpendicular to the plane of the cut are allowed. 

— * 

In these cases, / is again identified with the relative displacement of the centers 
of the two halves of the crystal. When symmetry constraints are absent, atomic 
relaxations in all directions are allowed. In these cases, the correct value of the 
relative displacement / is obtained from the relaxed in-plane coordinates of atoms 
on the two planes immediately adjacent to the cut, rather than from the relative 
displacement of the centers of the two halves of the crystal. Therefore, for a given 

— * 

relative displacement fc between the centers of the two halves of the crystal, atomic 

relaxation not only reduces the energy but also defines the actual displacement / 
(generally / 7^ fc), which is the value relevant to the GSF energy surface. 
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3 GSF Energies 

We first consider tlie (111) plane wliicli is the natural cleavage plane of silicon. For 
this plane of the diamond lattice, there are two distinct ways to cut the crystal: 
the shuffle set and the glide set. The difference between these two ways is shown 
in Fig. 1(a). For the shuffle set, the vertical distance between the two adjacent 
atomic planes immediately above and below the cut is equal to the bond length; 
only one bond per atom on either side of the cut is broken during the block shearing 
process. The corresponding interplanar distance for the glide set is only 1/3 of the 
bond length; three bonds per atom on either side of the cut are broken in the slip 
along this plane. The calculated GSF energy surfaces for the corresponding cuts 
in rigid block shearing with all the atoms in the two respective half crystals held 
fixed during shearing, are shown in Fig. 2(a) and 3(a) respectively. The effects 
of atomic relaxation on the GSF energy surfaces are shown in Fig. 2(b) and 3(b) 
for the shuffle and glide sets. Because the distance between the two atomic planes 
adjacent to the cut is sufficiently large in the shuffle set, the relaxation effects on 
this energy surface are very small and hardly noticeable on the scale of jus, which 
is of order 2 eV [compare Fig. 2(a) and 2(b)]. 

The effect of atomic relaxation is much more pronounced in the glide set. In 
Fig. 1 we show the atomic positions on a (110) plane, during the shearing on the 
glide plane along the < 121 > direction, for unrelaxed structures. Fig. |I|(a) is the 
original ideal structure. Fig. ^^b) is the structure for the unstable stacking fault 
configuration. Fig. |I](c) is the stable stacking fault configuration for the glide set. 
There is essentially no difference between the stable stacking fault structure and 
the ideal structure, as far as the coordination numbers of the two sets of atoms 
adjacent to the cut are concerned. For this reason, the energy of the stable stacking 
fault configuration is very low (see Fig. 3). For displacements beyond the stable 
stacking fault configuration along the < 121 > direction, the atoms belonging 
to the same color (white or black) immediately adjacent to the cut start getting 
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very close to each other as can be seen from Fig. |l](d), (e) and (f). Relaxation 
is restricted by symmetry to the direction perpendicular to the plane of the cut 
for configurations [3(b), (c) and (e). Especially in Fig. |l|(e), the atoms are within 
a very small distance, almost on top of each other. This kind of arrangement of 
atoms costs a significant amount of energy due to the strong repulsive interaction 
between the ionic cores. Configuration |l](e) actually corresponds to the highest 
energy for the glide cut. The relaxed energy surface on the glide set has the same 
energy scale as that for the shuffle set. The reason for this dramatic reduction in 
energy [compare Fig. 3(a) and 3(b)] is that through relaxation the atoms can avoid 
coming very close to each other during the sliding process. 

The results for the (100) plane before and after atomic relaxation of the atomic 
coordinates are shown in Fig. 4(a) and 4(b) respectively. The origin for the high 
energy barrier along the <011> path is similar to the situation for the glide set, 
where atoms adjacent to the cut come very close to each other. Therefore, the 
relaxation of the atomic coordinates is expected to reduce the energy barrier along 
this direction significantly. As is shown in Fig. 4, the topology of the energy surface 
is different after atoms are allowed to relax. One point deserves further attention: 
the <011> path, which is the path containing the highest energy barrier before 
relaxation, has become the energetically favorable path for the gliding process on 
the (100) plane after taking into account the effects of atomic relaxation. 

The calculated values of 'jus both before and after the relaxation of atomic 
coordinates for the (111) and (100) planes of silicon are summarized in Table I. 
The expressions that give the relaxed GSF energy surfaces, fitted by sinusoidal 
expansions, are given in the Appendix. 

Based on Rice's theory, a direct comparison of the values of 'jus indicates that 
it is energetically more favorable to nucleate dislocations on the shuffle set of the 
(111) plane under zero temperature and zero pressure conditions. The changes in 
the energy surface due to atomic relaxation will also affect the estimates of the 
free energy. Following Kaxiras and Duesbery (1993), we define the free energy per 
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unit area associated with a particular slip process as: 

F = lus-T^-P^. (1) 

where 5* is the entropy and /W the volume relaxation at the saddle point. AV is 
obtained by minimizing the energy of the saddle point configuration with respect 
to the magnitude of the slab thickness. We use Vineyard's transition state theory 
(Vineyard (1957)) as outlined by Kaxiras and Duesbery (1993), to estimate the 
entropy S for the relaxed energy surface. The condition that the preferred slip 
plane changes from shuffle to glide is then given by F^'*'^^ = p^^^uffie^ T]ie {P,T) 
values that satisfy this condition are shown in Fig. |^. As anticipated in the work 
of Kaxiras and Duesbery (1993), the larger reduction in the energy associated 
with relaxation of the glide cut results in an overall shift of the F^^^'^^ = F^^^^^'-'^ 
line towards higher pressure and lower temperature, compared to the unrelaxed 
results. For example, the stress needed to change from the shuffle to the glide set 
at room temperature (300° K) corresponds to an external pressure of only 15 kbar 
as calculated from the results with the relaxed energy surface, while a pressure as 
high as 53 kbar would be predicted if the unrelaxed energy surface were used. 

However, even the results obtained from the relaxed energy surface calcula- 
tions should not be taken literally for the following reasons: (1) Other factors 
(such as the electronic degrees of freedom), which are neglected in the present en- 
tropy calculation, could make a significant contribution to the free energy. (2) The 
comparison of free energies is based on Rice's theory in which the energy associ- 
ated with the creation of surface during dislocation emission is neglected; further 
investigations are needed to get an estimate of such effects (see for instance Xu, 
Argon and Ortiz (1995)). A detailed discussion on surface effects in dislocation 
nucleation for silicon will be published elsewhere (Juan, Kaxiras, and Sun (1995)). 
(3) Effects related to reconstruction of the dislocation core are also neglected in this 
theoretical framework. For a stiff covalent material, core reconstruction effects are 
very important (see Bulatov, Yip and Argon (1995)), and are expected to change 
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the values of '-jus- All these approximations may affect the exact numbers for the 
transition from shuffle to glide dominance, even though we expect the qualitative 
picture to remain the same. 



4 Peierls-Nabarro Model for Planar Dislocation 

By combining our results for the GSF energy surfaces with the Peierls-Nabarro 
model (Peierls (1940), Nabarro (1947)), we can further analyze the properties of 
dislocations on these planes within continuum elastic theory. We give here a brief 
description of the Peierls-Nabarro model which is based on the premise that a 
dislocation can be thought of as a continuous distribution of infinitesimal crystal 
misfits on the glide plane (Kroupa and Lejcek (1972), Hirth and Lothe (1982), Joos 
et al. (1994)). In the following we treat the distortion due to the dislocation as a 
scalar function of position f{x) for simplicity, although the true vectorial character 

— * 

of / should be considered in a general treatment. For each point on the ghde plane 

at a distance x' from the dislocation line, there is a corresponding infinitesimal 
Burgers vector df = df /dx\x'dx' , where df/dx is defined to be the dislocation 
density p{x) and f{x) is the disregistry at the point x, satisfying the condition: 



with b the Burgers vector of the dislocation. Prom the elastic model, the stress 
due to a dislocation along the direction of the Burgers vector on the slip plane is 

given by 

(^bn = 7^ , (<3j 

where x is the distance between the dislocation and the point where the stress is 
evaluated, and n is the normal to the glide plane. The constant K depends on the 

clastic properties of the crystal and the dislocation type. Its value for an isotropic 
solid is given by Joos et al. (1994): 





sm\e) 
(1^ 



+ cos'(^) , 



(4) 
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where fi and u are the shear modulus and Poisson's ratio respectively, and 6 is the 
angle between the dislocation line and the Burgers vector. Therefore, the stress 
due to the existence of the dislocation for a point x on the glide plane is: 

I '-dx'. (5) 



r^-^dx' = F,{f{x)). (6) 

J -co X — X' 



-oo X — X' 

Meanwhile, there is another stress which is due to the local disregistry of the crystal 
and tends to restore the crystal. This crystal restoring stress is a periodic function 
of the displacement and can be represented as Fb(/). Equilibrium is attained when 
the stresses due to the two terms are in balance with each other as expressed in 
the following equation, known as the Peierls-Nabarro equation, 

K r+°° p{x') 

2tT J-oo X — x' 

The dislocation profile f{x) and the dislocation density p{x) — df /dx can be deter- 
mined by solving this equation with the normalization condition p{x)dx — b. 
In the original Peierls-Nabarro model, a simple sinusoidal form is assumed for the 
restoring stress, 

W(^))=i^^a.sin(?^), (7) 
with F^ax the maximum stress. This assumption leads to the analytic solution 

/(^) = -tan-i7 + ^- (8) 

The parameter ( — Kb/AnFjnax can be viewed as the width of the dislocation 
within this model, since the value of the dislocation density p at this point is exactly 
one half its value at a; = 0. Within the GSF approach the crystal restoring force 
can be obtained by simply taking the derivative of the GSF energy surface 
7(/) with respect to the lattice distortion / (Kroupa and Lejcek (1972), Vitek 
and Yamaguchi (1973), Yamaguchi and Vitek (1975)). For a general functional 
form of F, there is no analytic solution for f{x). We follow the procedure of Joos 
et al. (1994) to solve the Peierls-Nabarro equation numerically. This is done by 
expressing the disregistry vector f{x) as a series 

fix)^^±a,t^n-'^+^-, (9) 
^ i=i Q 2 
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with the parameters a^, Xi to be determined, subject to the normahzation 
condition 

E«^ = i, (10) 

1=1 

since the total Burgers vector should be equal to b. By substituting the above 
displacement formula into the two sides of the Peierls-Nabarro equation, we get 
the forces F{ai,(i,Xi,x) and F'{ai,(i,Xi,x) corresponding to the contributions 
from dislocation distribution and crystal distortion respectively. The difference 
between these two forces \F{ai,Q then minimized by 

varying the parameters a^, Q and Xj so that the Peierls-Nabarro equation is satisfied 
numerically. We find that a numerical solution is feasible by retaining only three 
terms in Eq. (||). 

We display in Fig. ^ the GSF energies and forces obtained from our first- 
principles calculations with the full relaxation of the atomic coordinates, along the 
directions which are relevant to the dislocations we will consider. For the case of 
the < 121 > direction on the glide plane, only the portion ranging from the ideal 
structure to the stacking fault configuration is shown, since a partial dislocation 
will be formed along this direction. One feature of these curves deserves further 
comment: not all the restoring forces can be well described by a simple sinusoidal 
from (see also the recent work of Xu, Argon and Ortiz (1995)). Especially in the 
case of the (100) plane, the force curve is very flat when the displacement is approx- 
imately one half of the full Burgers vector. The significance of this deviation will 
be discussed in detail later, when we consider the profiles of various dislocations. 

The disregistry vector for full dislocations on the shuffle set of the (111) plane, 
the glide set of the (111) plane, and the (100) plane are displayed in Fig. |^ (a), 
(b), and (c) respectively. We note that the dislocations with Burgers vector along 
the < 101 > direction on the glide plane (glide-60° and glide-screw) are more 
concentrated compared to the other dislocations on the (111) plane (this effect will 
become more obvious when we compare the dislocation density). This difference 
in the distribution density can be understood by considering the magnitude of 
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the crystal restoring force for these dislocations. A qualitative estimate can be 
obtained by considering the half width of the dislocation ( within the classical 
model, which is inversely proportional to the magnitude of the crystal restoring 
force. This force is largest in the < 101 > direction of the glide plane, as is evident 
by comparing Figures 2, 3 and 4. 

In Fig. D we show representative dislocation densities for the planes we have 
considered, including (a) the 60° dislocation on the shuffle set, (b) the 60" dis- 
location on the glide set, and (c) the 90° dislocation on the (100) plane. We 
have displayed the results obtained from both the relaxed (upper panel) and the 
unrelaxed GSF energy surfaces (lower panel) to examine the effects due to the 
relaxation of the atomic coordinates. For comparison we have also displayed the 
results from the classical solution (dashed curves), obtained by assuming a sinu- 
soidal form for the crystal restoring force with its maximum equal to the value 
obtained from the first-principles calculations. For the case of the 60" dislocation 
on the shuffle set, it is apparent that the use of the relaxed GSF energy surface 
does not cause significant change on the dislocation profile both qualitatively and 
quantitatively. This is a consequence of the fact that atomic relaxation does not 
change the GSF energies along this direction in any significant way. For the 60" 
dislocation on the glide set, the use of the relaxed energy surface makes the dis- 
location profile considerably smoother and wider but the general shape remains 
unchanged. This is due to the fact that the relaxation reduces the magnitude of 
the crystal restoring force but does not change its functional form significantly. For 
the case of the 90" dislocation on the (100) plane, the dislocation profile is changed 
qualitatively with the use of the relaxed energy surface. This reflects the fact that 
the crystal restoring force for this plane is significantly different from the sinusoidal 
function as we mentioned above. The double peak indicates the dissociation of the 
full dislocation into two fractional dislocations (Vitek and Kroupa (1969)). A brief 
summary of the obtained dislocation properties is given in Table II. 

Another useful quantity to consider is the energy barrier associated with dislo- 
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cation motion. Within the Peierls-Nabarro model, this energy barrier is called the 
Peierls energy and is defined as the amplitude of the variation of the misfit energy 
on the glide plane as the position of the dislocation line moves. With the obtained 
dislocation profiles, we can calculate the misfit energy across the glide plane as a 
function of the position of the dislocation line u following the definition given by 
Vitek and Yamaguchi (1973) and Joos et al. (1994). To be consistent with the 
Peierls-Nabarro model, where the discreteness of the crystal is taken into account 
at the glide plane, the misfit energy W{u) is defined as the sum of all the misfit 
energies between pairs of atom rows, obtained from the GSF energy at the local 
disregistry: 

+00 

Wiu)= J2 lifima'-u))a', (11) 

m=—ao 

where a' is the distance between adjacent atomic planes in the direction perpen- 
dicular to the dislocation line. This expression meets two important requirements: 
First, it has the correct period a' of planes in the crystal, 

W{u + a') =W{u). (12) 

Second, in the limit of a very narrow dislocation, it reproduces the correct maxi- 
mum 7usa'. The other quantity which might be of interest is W{a' /2) where the 
minimum of the misfit energy function occurs. Since this quantity measures the 
nonelastic part of the energy of the dislocation, it provides a qualitative estimate 
of the core energy. The stress associated with the energy function W{u) can then 
be defined as: 

, , 1 dW(u') , , , 

The maximum of this stress function is the Peierls stress, needed to move the 
dislocation. The values of the energy function W{u) and the stress function (t{u) 
for the (111) plane are shown in Fig. ^ (a) and (b). It is obvious that the dislocation 
with Burgers vector along the < lOl > direction on the glide set (labeled Glide-60° 
in Fig. y (b)) has larger Peierls energy and Peierls stress. This is expected since 
we know from our GSF energy surface calculations that the energy barrier along 
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that direction is higher than those along the other two vectors. Fig. ^(c) is the 
energy and stress function for the case of the (100) plane. The different behavior 
between the 90*^ and the screw dislocations on the (100) plane is primarily due to 
the difference in the values of the constant K. The Peierls energy, stress, and core 
energy of various dislocations are summarized in Table III. For the (111) plane, it 
appears that the dislocations belonging to the shuffle set are the easiest to move 
under zero temperature and zero pressure conditions. However, as we discussed 
earlier, including the effects of temperature, pressure, surface creation and core 
reconstruction can change the picture. Finally, we would like to mention that 
there is experimental evidence that on the glide set the 90° partial is more mobile 
than the 30° partial (Wessel and Alexander (1977), Alexander, Gottschalk, and 
Kisielowski-Kemmerich (1985), Grosbras, Demenet, Garem, and Desoyer (1984), 
Demenet, Grosbras, Garem, and Desoyer (1989)), which is consistent with the 
results of our calculations, namely a{u) is lower for the glide-90°(p) than the glide- 
30° (p) dislocation, see Fig. 9(b). 



5 Summary 

In summary, we have performed first-principles calculations to obtain the GSF 
energy surfaces for both the shuffle and glide sets of the (111) plane, as well as the 
(100) plane of silicon. We showed that for the glide set and the (100) plane, the 
effects on the GSF energy surfaces due to the relaxation of the atomic coordinates 
are significant. The unstable stacking fault energies '-^us for these planes were de- 
termined from our calculations. By combining these values with Rice's criterion for 
dislocation nucleation, the shuffle set appears to be favored for dislocation emis- 
sion under zero temperature and zero pressure conditions. A qualitative account 
of entropy effects was attempted, based on Vineyard's transition state theory. By 
comparing the free energies of the shuffle and glide sets, we find that either set can 
dominate under different thermodynamical conditions. Quantitative free-energy 
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comparisons should also take into account additional entropy effects, and the ef- 
fects of surface creation and core reconstruction. 

The GSF energy surfaces were then combined with the Peierls-Nabarro model 
to investigate the dislocation properties of silicon. The crystal restoring forces are 
obtained directly from our first-principles GSF energy surfaces. We demonstrated 
the importance of using the relaxed energy surface for the calculation of these 
properties. The differences between dislocation profiles obtained with relaxed GSF 
and unrelaxed GSF energy surface are significant in certain cases. The Peierls 
energy and Peierls stress were also calculated within the framework of the Peierls- 
Nabarro model. 

We would like to conclude with a brief discussion on the approximations in- 
volved in the Peierls-Nabarro model: (1) The dislocation is assumed to be planar 
within this model. We expect that the errors introduced with this approximation 
for a stiff material hke sihcon should be minimal. (2) The response from the crystal 
is treated within elastic theory, which may not apply for a very narrow core situa- 
tion, where non-elastic effects are expected to be important. (3) The Peierls stress 
concept assumes that the only relevant quantity in determining the stress is the lo- 
cal disregistry, which is true only when the dislocation density is very smooth. (4) 
The dislocation line is assumed to move as a rigid object, which is unlikely (see e.g. 
the recent work of Bulatov, Yip and Argon (1995) on this issue). Accordingly, we 
suggest that these results should not be taken literally. Rather, our investigation 
provides a qualitative comparative discussion of dislocation properties in silicon, 
that can serve as guide to more detailed studies. 
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Appendix 



We used sinusoidal expansions to fit the calculated GSF energy surfaces in order to 
facilitate the computation of dislocation properties. The expansions were chosen 
so that they satisfy the underlying translational symmetries of the lattice. We 
have checked the numerical values and found that the numbers obtained from the 
fittng formula reflect the underlying rotational symmetry with sufficient accuracy. 
Specifically, we use the following expression: 

/ \ ^ .2?™ , ,27rm , „ ,27m , . ,27rm , 
7(a;, y) = ^ A^m cos( x) cos( y) + cos( x) sm( y) 

n,m '^1 "^2 Ol O2 

^ . ,27rn , . ,27rm , ^, 
+Cnm sm( x) sm( y) (14) 

where oi and 02 are the repeat distances in the (x, y) plane. In terms of the lattice 



constant oq of bulk Si, these are given by a\ — oq/v^, 02 = Oq x ^3/2 for the (111) 
plane cuts (both shuffle and glide), and oi = 02 = oo/v^ for the (100) plane. The 
{x^y) directions in the (111) plane cuts correspond to the < 101 > and < 121 > 
crystallographic directions, whereas in the (100) plane cut they correspond to the 
< Oil > and < Oil > crystallographic directions. The coefficients of the terms 
retained in the above expansion are given in Table IV. 
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Table I 





No Relaxation 


Atomic Relaxation 
at ideal volume 




location 


'Jus (J/m2) 


location 


(J/m2) 


(lll)-Shuffle 


am 


1.84 


i[ioi] 


1.81 


(lll)-Glide 




2.51 


^[121] 


2.02 


(100) 




2.97 


i[oii] 


2.15 



TABLE I: The unstable stacking fault energy, 7„s, obtained for the shuffle set 
and glide set of the (111) plane, as well as the (100) plane of silicon. The location 
indicates the position of relative displacement / where the jus occurs. 



Table II 



dislocation 


Shuffle 


Glide 


Glide partial 


(100) 


60° 


Screw 


60° 


Screw 


30° 


90° 


90° 


Screw 


X(10ii^^) 

V cm'' ' 


8.02 


6.37 


8.02 


6.37 


6.92 


8.58 


9.04 


6.375 


'''maa; (10 ) 


1.49 


1.49 


4.29 


4.29 


2.78 


2.78 


2.06 


2.06 


C(A) 


1.08 


0.92 


0.46 


0.37 


0.77 


0.92 


2.15 


1.54 


b(A) 


3.84 


3.84 


3.84 


3.84 


2.22 


2.22 


3.84 


3.84 



Table II: Quantities related to the properties of dislocations for the shuffle and 
glide sets on the (111) plane, as well as the (100) plane of silicon. The meaning of 
K is given in Eq. (^), T^ax is the maximum value of d-^us^D/df for the directions 
we are interested in, C, is the half width of the dislocation, and h is the Burgers 
vector. 
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Table III 



dislocation 


Shuffle 


Glide 


Glide partial 


(100) 


60° 


Screw 


60° 


Screw 


30° 


90° 


90° 


Screw 


WpieVA) 


0.149 


0.183 


0.842 


0.898 


0.287 


0.246 


0.112 


0.296 




0.046 


0.062 


0.399 


0.504 


0.176 


0.139 


0.032 


0.081 


W{a'/2) 


0.467 


0.338 


0.345 


0.238 


0.174 


0.246 


0.549 


0.274 



Table III: The Peierls energy Wp , Peierls stress Up, and the estimated core 
energy W{a'/2) obtained from the Peierls-Nabarro model for various dislocations 
on the (111) and (100) planes of silicon. 

Table IV 





(lll)-shuffle 


(lll)-glide 


(100) 


(n,m) 


A 






A 




c 


A 


(0,0) 


0.9789 






2.2683 






1.8463 


(0,1) 














-0.6728 


(0,2) 


-0.3857 




0.0625 


-0.7059 




0.7024 




(0,4) 


0.0180 




-0.0200 


0.2570 




-0.3473 




(1,0) 














-0.8443 


(1,1) 


-0.7714 


-0.1231 




-1.3609 


-1.8219 




-0.2670 


(1,2) 














-0.1493 


(1,3) 


0.0950 


-0.0010 




-0.6291 


0.1218 






(2,0) 


0.0470 






-0.0126 






0.1471 


(2,1) 














0.0087 


(2,2) 


0.0350 


0.0390 




0.4412 


0.0890 




-0.0485 


(2,4) 








-0.2579 


0.2698 






(3.1) 


-0.0180 


-0.0070 













Table IV: The non-zero coefficients for the expansion of the 7 surface in terms 
of sinusoidal functions for the shuffie and glide (111) cuts, and for the (100) cut. 
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Figure 1: The atomic structures involved in the shearing process, in which one 
half of the diamond crystal is sheared rigidly with respect to the other half along 
the < 121 > direction on the glide set of the (111) plane of silicon. Both the 
glide (G) and shuffle (S) cuts are indicated in part (a) by dashed lines. The two 
different shadings of the atoms represent two different (110) atomic planes, (a) 
the ideal structure, (b) the unstable stacking fault configuration, (c) the stacking 
fault configuration, (d), (e), and (f) are intermediate geometries, when the upper 
half is sheared further along the < 121 > direction, (e) is the geometry with the 
highest energy. 

Figure 2: The GSF energy surface: (a) before and (b) after relaxation, for the 
shuffle set of the (111) plane. 

Figure 3: The GSF energy surface: (a) before and (b) after relaxation, for the 
glide set of the (111) plane. 

Figure 4: The GSF energy surface: (a) before and (b) after relaxation, for the 
(100) plane. 

Figure 5: The shuffle set / glide set transition phase diagram on the (P, T) plane. 
Results including atomic relaxation (solid line) and before relaxation (dashed line) 
are shown for comparison. For (P, T) values below (above) the hne, the ghde 
(shuffle) has lower unstable stacking free energy. 

Figure 6: The fully relaxed GSF energy curve and the corresponding stress for the 
relevant directions on the shuffle, glide and (100) planes of silicon. 

Figure 7: (a) The disregistry vector obtained for dislocations on the shuffle set of 
the (111) plane; (b) the disregistry vector for dislocations on the glide set of the 
(111) plane; (c) the disregistry vector for dislocations on the (100) plane. 
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Figure 8: The dislocation density obtained from our calculations for (a) the 60° 
dislocation on the shuffle set, (b) the 60^^ dislocation on the ghde set, and (c) the 
90° dislocation on the (100) plane. Results from both relaxed (upper panel) and 
unrelaxed (lower panel) GSF energy surfaces are displayed for comparison. The 
dashed curves are the classical solutions (see text). 



Figure 9: (a) The energy function W{u) with respect to the displacement vector 
u (upper panel) and the corresponding stress a{u) as a function of u (lower panel) 

for dislocations belonging to shuffle set. (b) Same as in (a), for dislocations on the 
glide set. (c) Same as in (a) for dislocations on the (100) plane. 



23 



